Note
Click here to download the full example code
Drill-hole spacing from virtual grids¶
A mine plans quarterly production from 40 × 40 m areas of Walker Lake and needs each quarter's grade within ±15 % at 90 % confidence. Which drill spacing delivers that? Wilde (2010) answers it by drilling simulated deposits: take realizations conditioned to today's data as the truth, sample them on regular grids, simulate again from those samples alone, and measure the uncertainty of each production volume. You run his virtual grids at 5 to 20 m, read the spacing each criterion needs, and check the result against a reference study of the same data.
import sys
import warnings
from pathlib import Path
HERE = Path(__file__).parent if "__file__" in globals() else Path.cwd()
sys.path.insert(0, str(HERE.parents[1]))
warnings.filterwarnings("ignore", "Samples sharing")
import boitata as bt
import matplotlib.pyplot as plt
import numpy as np
from common import ACCENT, GRAY, HIGHLIGHT, INK, LIGHT, save
What you'll learn
- How Wilde's virtual grids turn a simulation model into an uncertainty-versus-spacing curve.
- What the maximum expected error (MEE) of a production volume measures, and how to read the share of volumes that meet ±15 % at 90 %.
- Why the required spacing depends on the volume: a quarter needs much denser drilling than a year.
Prerequisites: simulation at block support
for window= and SimulationSummary.relative_error.
The data¶
The 470 samples of V (ppm), declustered with cell weights, give the normal-score variogram: a nugget of 0.37 and
one spherical structure of range 103 m along azimuth 170°, a third of that across. Turning bands with 500 bands
simulate on 2 m nodes. The production volumes are moving windows. For a quarter, blocks= averages the nodes to
8 m panels and window=(40, 40) gives each panel the mean of the 5 × 5 panels around it, a 40 × 40 m window;
for a year, 16 m panels and window=(80, 80). A window spans an odd number of panels, so the panel side is a
fifth of the window. The study counts only windows that lie wholly inside the area.
samples = bt.datasets.walker_lake()
xy, v = samples.coords, samples["V"]
weights = bt.cell_declustering(samples, "V", sizes=np.arange(2.5, 102.5, 2.5)).weights
scores = bt.NormalScore().fit(v, weights=weights).transform(v)
azimuths = (170, 260)
directional = [bt.experimental_variogram(xy, scores, 10, 120, azimuth=a) for a in azimuths]
fitted = bt.Variogram.fit_directional(directional, [(a, 0) for a in azimuths], rotation=[170, 0, 0])
gaussian = fitted.standardized()
print(gaussian)
nodes = bt.BlockModel(origin=(1, 1), size=(2, 2), count=(130, 150))
panels = {
side: bt.BlockModel(origin=(1, 1), size=(side / 5,) * 2, count=(1280 // side, 1480 // side))
for side in (40, 80)
}
tb = bt.TurningBands(gaussian, bands=500, search=bt.Search(radius=100, max_samples=24)).fit(
xy, v, weights=weights
)
def inside(side):
"""Panels whose window of `side` meters lies inside the panels."""
c = panels[side].coords[:, :2] - 1
extent = np.array([1280 // side, 1480 // side]) * side / 5
return ((c >= side / 2) & (c <= extent - side / 2)).all(axis=1)
print(f"{len(nodes.coords)} nodes, {inside(40).sum()} quarterly and {inside(80).sum()} yearly windows inside")
Out:
Variogram(nugget=0.3712292816104307, structures=[Structure("spherical", sill=0.6287707183895693, range=102.93174851181992)], rotation=(170.0, 0.0, 0.0), ratios=(0.34469006793164536, 1.0))
19500 nodes, 924 quarterly and 168 yearly windows inside
Step 1: Pick the truths
Each truth is a realization conditioned to the 470 samples, so it honors what is known today and varies
where nothing is known. Twenty candidates are simulated and averaged to the quarterly windows;
select_realizations keeps three that span them, the medoids of a k-medoids clustering of their window
grades. A truth index is the realization's number in simulate(seed=101), which spacing_study simulates
again.
candidates = tb.simulate(nodes, n=20, seed=101, blocks=panels[40], window=(40, 40), keep=True)
truths = bt.select_realizations(candidates, 3, seed=0)
means = candidates.realizations[:, inside(40)].mean(axis=1)
print(
f"candidate means {means.min():.0f} to {means.max():.0f} ppm; truths {truths.tolist()}: "
+ ", ".join(f"{means[k]:.0f}" for k in truths)
+ " ppm"
)
Out:
candidate means 302 to 335 ppm; truths [9, 19, 5]: 318, 307, 334 ppm
Step 2: Drill each truth on virtual grids
spacing_study lays vertical holes on a square grid over the nodes for each spacing, reads each truth at
the holes, refits a copy of the simulator on those samples alone and simulates 50 realizations of the
windows. The model is 2D, so composite_length=np.inf gives one sample per hole. Each table row is one
window of one truth at one spacing.
spacings = [5, 7.5, 10, 12.5, 15, 20]
study = {
side: bt.spacing_study(
tb,
nodes,
spacings=spacings,
truths=truths,
seed=101,
n=50,
blocks=panels[side],
window=(side, side),
composite_length=np.inf,
)
for side in (40, 80)
}
for side, table in study.items():
table = table.filter(inside(side)[np.asarray(table["row"], dtype=int)])
study[side] = table
print(f"{side} m windows: {table.num_rows} rows, columns {', '.join(table.column_names)}")
Out:
40 m windows: 16632 rows, columns plan, spacing, realization, row, truth, mean, mee, cv, precision_0.15, error, covered
80 m windows: 3024 rows, columns plan, spacing, realization, row, truth, mean, mee, cv, precision_0.15, error, covered
The math
For a window with realizations \(z^{(1)}, \dots, z^{(n)}\), mean \(\bar z\) and 5th and 95th percentiles \(q_{05}, q_{95}\), the maximum expected error (Koppe et al., 2017) is [ \mathrm{MEE} = \frac{q_{95} - q_{05}}{2\,\bar z}. ] A window with MEE ≤ 0.15 has 90 % of its realizations within about ±15 % of the mean.
Step 3: MEE per spacing
The box plot shows the spread of MEE over the quarterly windows of the three truths; the right panel counts the windows that meet ±15 % at 90 % for both volumes.
spacing_of = {side: table["spacing"] for side, table in study.items()}
share = {
side: [np.mean(t["mee"][spacing_of[side] == s] <= 0.15) for s in spacings] for side, t in study.items()
}
print("spacing median MEE 40 m share ≤ 15 %: 40 m 80 m")
for i, s in enumerate(spacings):
median = np.median(study[40]["mee"][spacing_of[40] == s])
print(f"{s:7g} {median:13.1%} {share[40][i]:17.1%} {share[80][i]:5.1%}")
fig, (a, b) = plt.subplots(1, 2, figsize=(11, 4.2), layout="constrained", width_ratios=[1.2, 1])
boxes = a.boxplot(
[100 * study[40]["mee"][spacing_of[40] == s] for s in spacings],
positions=spacings,
widths=1.4,
whis=(10, 90),
showfliers=False,
patch_artist=True,
medianprops={"color": INK, "lw": 1.5},
)
for box in boxes["boxes"]:
box.set(facecolor=LIGHT, edgecolor=GRAY)
a.axhline(15, color=HIGHLIGHT, ls="--", lw=1.2)
a.set(
xlabel="Virtual grid spacing (m)",
ylabel="MEE of the window (%)",
xticks=spacings,
xticklabels=[f"{s:g}" for s in spacings],
title="40 × 40 m windows: MEE (P10–P90)",
)
for side, color in ((40, ACCENT), (80, HIGHLIGHT)):
b.plot(spacings, 100 * np.array(share[side]), "o-", color=color, label=f"{side} × {side} m")
b.axhline(90, color=GRAY, ls="--", lw=1)
b.set(
xlabel="Virtual grid spacing (m)",
ylabel="Windows with MEE ≤ 15 % (%)",
ylim=(0, 102),
title="Share meeting ±15 % at 90 %",
)
b.legend()
save(fig, "virtual_grids")
Out:
spacing median MEE 40 m share ≤ 15 %: 40 m 80 m
5 10.2% 89.1% 100.0%
7.5 14.6% 53.1% 99.8%
10 18.4% 22.4% 94.2%
12.5 22.1% 9.8% 81.3%
15 26.5% 3.8% 48.0%
20 32.1% 0.0% 14.3%
The median quarterly MEE grows from 10.2 % at 5 m to 32.1 % at 20 m. At 5 m, 89.1 % of the quarterly windows meet the criterion; at 7.5 m, 53.1 %; at 10 m, 22.4 %. A yearly window covers four times the area, and 94.2 % of them still meet it at 10 m.
Step 4: Required spacing
uncertainty_curve bins the windows by spacing and gives the P50 and P90 of MEE in each bin;
required_spacing interpolates where a quantile crosses 15 %. The P90 crossing is the spacing at which 90 %
of the windows meet ±15 %; the P50 crossing, at which the median window does.
edges = [4, 6, 8.5, 11, 13.5, 17, 22]
curves = {side: bt.uncertainty_curve("spacing", "mee", bins=edges, data=t) for side, t in study.items()}
required = {}
for side, curve in curves.items():
with warnings.catch_warnings():
warnings.simplefilter("ignore")
required[side] = {q: bt.required_spacing(curve, column=q) for q in ("P50", "P90")}
print(
f"{side} m windows: P90 MEE at 5 m {curve['P90'][0]:.1%}; required spacing for the median window "
f"{required[side]['P50']:.1f} m, for 90 % of windows {required[side]['P90']:.1f} m"
)
fig, axes = plt.subplots(1, 2, figsize=(11, 3.8), layout="constrained")
for ax, (side, curve) in zip(axes, curves.items(), strict=True):
bt.plot.uncertainty_curve(curve, required=required[side]["P90"], ax=ax)
ax.set(title=f"{side} × {side} m windows", xlabel="Virtual grid spacing (m)", ylabel="MEE")
save(fig, "required_spacing")
Out:
40 m windows: P90 MEE at 5 m 15.3%; required spacing for the median window 7.7 m, for 90 % of windows nan m
80 m windows: P90 MEE at 5 m 7.1%; required spacing for the median window 14.8 m, for 90 % of windows 11.4 m
The yearly volume needs holes every 11.4 m for 90 % of its windows, and its median window meets ±15 % up to
14.8 m. The quarterly volume misses the 90 % mark on every grid tried: its P90 MEE is 15.3 % at 5 m, so
required_spacing returns NaN. Its median window crosses 15 % at 7.7 m.
Key idea
The required spacing belongs to a volume and a confidence, not to the deposit. Averaging over a larger volume cancels more of the local error, so a year tolerates holes twice as far apart as a quarter.
Check before you move on
First, covered says whether the truth of a window lies within its MEE of the simulated mean; a calibrated
model covers about 90 % of the windows. Second, a reference study of the same data (Wilde's method on 2 m
nodes, three truths, 50 realizations) found 83.6 % of the quarterly windows within ±15 % at 5 m and 24.2 % at
10 m. Its holes sat half a spacing in from the edge, between the nodes. planned_drillholes with that
offset builds the same grids, and spacing_study(plans=) drills them on ten truths.
for side, table in study.items():
print(f"{side} m windows covered: {np.mean(table['covered']):.1%}")
plans = {f"{s:g}": bt.planned_drillholes(nodes, s, offset=(s / 2 - 2, s / 2 - 2)) for s in (5, 10)}
reference = bt.spacing_study(
tb,
nodes,
plans=plans,
truths=10,
seed=101,
n=50,
blocks=panels[40],
window=(40, 40),
composite_length=np.inf,
)
reference = reference.filter(inside(40)[np.asarray(reference["row"], dtype=int)])
plan = reference["plan"]
for name in plans:
meets = reference["mee"][plan == name] <= 0.15
truth_of = reference["realization"][plan == name]
per_truth = [np.mean(meets[truth_of == k]) for k in range(10)]
print(
f"holes between nodes, {name} m: {meets.mean():.1%} of quarterly windows; per truth "
f"{min(per_truth):.1%} to {max(per_truth):.1%}, standard deviation {np.std(per_truth, ddof=1):.1%}"
)
Out:
40 m windows covered: 86.8%
80 m windows covered: 86.0%
holes between nodes, 5 m: 84.5% of quarterly windows; per truth 78.4% to 90.6%, standard deviation 4.1%
holes between nodes, 10 m: 21.8% of quarterly windows; per truth 13.3% to 26.0%, standard deviation 3.7%
The MEE covers the truth in 86.8 % of the quarterly windows and 86.0 % of the yearly ones, a little short of 90 %. With the holes between nodes, 84.5 % of the quarterly windows meet ±15 % at 5 m and 21.8 % at 10 m, against the reference's 83.6 % and 24.2 %. The share changes from truth to truth: at 10 m one truth scores 13.3 %, another 26.0 %, a standard deviation of 3.7 points. The reference averaged three truths, so its share carries a standard error of about 2 points, and gaps of 0.9 and 2.4 points lie within Monte-Carlo noise. The grids of Step 2 put holes on node centers, and they score 89.1 % at 5 m: a node that holds a sample keeps its value in every realization, which narrows its window.
The decision¶
Drill on 11 m or closer for yearly planning: at 10 m, 94.2 % of the 80 × 80 m windows meet ±15 % at 90 %. A quarter needs holes closer than 5 m, which no resource drilling pays for; plan quarters from grade control. The learning curve reads the same relation from one simulation of the current data.
Total running time of the script: ( 0 minutes 22.637 seconds)
Download Python source code: example_14_07.py